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Abstract 


M.T.  Manry  and  J.K.  Aggarwal  recently  described  an  algorithm  for  use  in  the 

design  of  multi-dimensional  FIR  digital  filters  by  phase  correction.  As  they 

observe,  their  method  can  be  viewed  as  the  steepest  descent  method  for  minimizing 

a certain  function  f(x):  given  an  approximate  solution  x , a new  approximation 

n 

is  x ,,  - X +t  p where  p “-7f(x  ) and  t is  chosen  by  a simple  rule.  We 
n+1  n n*^n  *^n  n n j r 

derive  here  an  improved  rule  for  determining  t^  and  an  Improved  direction  p^ 

(essentially  the  Fletcher-Reeves  conjugate-gradient  direction).  The  resulting 
method  appears  to  be  two  to  three  times  as  fast  as  the  Manry-Aggarwal  method;  the 
additional  cost  is  primarily  in  storage,  which  roughly  doubles. 

Key  words:  FIR  digital  filters;  filter  design;  phase  correction; 

Manry-Aggarwal  method. 


1.  Introduction 

In  recent  papers  iManry  (1976),  Manry-Aggarwal  (1976)1,  M.T.  Manry  and  J.K. 
Aggarwal  have  proposed  a new  technique  for  the  design  of  multi-dimensional  FIR 
digital  filters.  The  reader  is  referred  to  those  papers  for  references,  applica- 
tions, and  comparisons  with  other  methods;  in  the  interest  of  brevity,  we  confine 
ourselves  here  to  presenting  dramatic  improvements  in  the  Manry-Aggarwal  method. 

For  simplicity  and  clarity,  we  follow  the  lead  of  [Manry  (1976) , Manry-Aggarwal 
(1976)1  by  presenting  our  discussion  in  terms  of  two-dimensional  filters;  general- 
ization is  obvious  and  straightforward. 


*Departments  of  Mathematics  and  of  Computer  Sciences  and  Center  for  Numerical 
Analysis  at  The  University  of  Texas  at  Austin.  Research  supported  in  part  by  the 
United  States  Office  of  Naval  Research  under  Contract  N00014-67-A-0126-0015; 
reproduction  in  whole  or  in  part  is  permitted  for  any  purposes  of  the  United  States 
government . 


2 


A two-dimensional  FIR  digital  filter  is  described  by  an  array  h of  filter 

spatial  coefficients  h for  o5m-M— 1 and  oin<N-l  for  integers  M,N. 

^ inn 

The  filter  produces  output  spatial  data  g^  from  input  spatial  data  d^^  according 
to 
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mn 


M-1  N-1 

''k£  ‘^m-k,n-£ 
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Corresponding  to  the  spatial  coefficient  array  h is  the  frequency-dong  in 
array  H of  coefficients  H for  o<m<M-l  and  o < n < N-1  defined  by  the 
discrete  Fourier  transform: 


(1.1) 
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where 

(1.3)  U - exp(-2nj/M) , V - exp(-2Tij /N) , j « >/  “1  . 


The  problem  addressed  by  Manry  and  Aggarwal  is  the  following.  It  is  desired 
to  design  a filter  whose  spatial  coefficients  h^  will  be  zero  except  for  a 
small  number  of  specified  values  of  m,n  and  whose  corresponding  frequency-domain- 
coefficient  amplitudes  [h^I  will  assume  (or  approximate)  prescribed  values 
for  o^k^M-1  and  oSH-N-1.  For  simplicity  of  presentation  it  is  assumed 
that  the  spatial  array  jj  is  truncated  in  the  sense  that  the  (possibly)  non-zero  spatial 
coefficients  h are  for  o5mSM,  -1,  o£n^N.  -1,  where  M.  iM,  N.  IN; 
generalization  is  obvious  and  straightforward.  In  this  case,  computation  of  the 
frequency-domain  array  H by  Equation  1.1  simplifies  to 
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Mj^-1  Nj^-1 
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2.  The  Manry-Aggarwal  method 

Manry  and  Aggarwal  first  present  a basic  iterative  step  for  improving  a 
spatial  array  to  an  array  which  comes  closer  to  having  the  desired 

frequency-domain-coefficient  amplitudes  Aj^.  Assuming  that  h^^  - o except  for 
o<m<Mj-l  and  o<n<Nj^-l,  we  compute  the  corresponding  frequency-domain  array 
by  Equations  1.1  or  1.4.  We  then  write  these  frequency-domain  coefficients 


(2.1) 


for  o<k<M-l,  o<«.<N-l. 


Recall  that  we  desire  to  have  for  oSk^M-l.  o^^^N-1.  forgiven 

^ a . 1 V 


f 1 

A^j^.  We  therefore  define  a new  frequency  domain  array  H with  the  same  phase 

as  for  but  with  the  correct  amplitude: 


(2.2) 


From  we  compute  the  corresponding  spatial  array  h^^^^  via  Equation  1.2. 

Since  we  cannot  in  general  expect  to  be  truncated,  that  is  to  satisfy  the 

requirement  that  e*oept  for  o^m^Mj-l  and  o^n^N^-l,  we  must 

truncate  to  obtain  h^  ^ : 


(2.3) 
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Ij^(l+*s)  o5m<M,-l,  o5n5N,-l 

rnn  I ^ 

o otherwise 


The  basic  iterative  step  then  takes  ••  h^^^^ . 

It  is  a simple  matter  to  see  that  h'  comes  closer,  in  a sense,  than 


M) 


to  achieving  the  desired  frequency-domain-coefficient  amplitudes  Aj^ . More 
iely,  given  any  truncated  spatial  i 
for  o5m5Mj^-l,  oSnlNj^-1,  we  define 


precisely,  given  any  truncated  spatial  array  h,  so  that  ° except  possibly 


(2.4) 


M-1  N-1 

f(h)  - E E (|H  I - V)' 
k-o  £-o 


where  H is  computed  from  h via  Equations  1.1  or  1.4.  Then  [Manry  (1976)J  we 
can  measure  the  improvement  in  Jj'  over  ^ by 


(2.5) 
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If  we  always  use  the  basic  iterative  step  » h^^^\  then  since  f is 

bounded  below  and  f decreases  at  each  step,  we  know  that  f(h^^^)  -f(h^^^^^) 
converges  to  zero;  by  Equation  2.5  this  implies  that  - h^^^  ((  converges  to 

zero,  and  hence  Joaniel  (1971^  the  sequence  either  converges  to  a unique 

* 

point  h or  it  has  a continuum  of  limit  points  (i.e.,  an  infinite  connected  set 

of  limit  points),  an  unlikely  occurence. 

In  Qtanry  (1976^  , another  important  interpretation  is  given  for  the  step  from 
( 1 ) ^ 

h to  h . The  function  f in  Equation  2.4  should  be  considered  a function 
of  the  M, N-  variables  h for  oSmiM, -1  and  o-n£N,-l.  It  can  then  be  shovm 
that  the  gradient  Vf  of  f is  given  by  (the  lo%#-order  components  of) 


(2.6) 


/ j \ A / 4 

Therefore  we  can  interpret  the  basic  iterative  step  from  J}'  ^ to  ^ as  a 


5 


step  of  length 


1 

2MN 


in  the  steepest  descent  direction 


for  f : 


(2.7)  . h^^^  + ^ 

Since  we  deduced  above  that  if  we  let  h^^  then  ||h^  h^  ^ ||  con- 

verges to  zero,  it  follows  from  Equation  2.6  that  Vf(h^^^)  converges  to  zero. 
Thus  we  can  analyze  the  convergence  of  the  basic  iterative  method. 


Proposition. 

Let  the  basic  iterative  method  be  used,  so  that  • h^^^\  Then  {h^ 

either  converges  to  a unique  point  or  has  a continuum  of  limit  points.  In  either 
case,  all  limits  h satisfy  7f(h)“o. 

Manry  and  Aggarwal  observed  that  in  practice  the  choice  ■ h^^^^  gave 

rather  slow  convergence;  this  simply  says  that  is  not  a particularly  good 

step  in  the  steepest  descent  direction.  They  therefore  propose  letting 


(2.8) 


,(1M)  . ^(i)  * ,S(1+S)  . h<l), 
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where  t^^  is  not  necessarily  chosen  as  * Their  choice  of  t^  attempts  to 

accelerate  convergence.  Specifically,  they  let  Cj  minimize  !lh^  ^ “ Ql 
c(h^^”^^  - 11^  with  respect  to  c,  giving  a simple  formula  for  c^;  they 

then  let  tj“l+c^  if  this  decreases  f,  but  t^^  ■ 2^  otherwise  or  if  i<2. 

Numerical  experiments  show  that  this  choice  gives  a dramatic  Improvement,  reducing 
f to  a given  value  in  half  the  work  as  for  the  basic,  unaccelerated,  method.  This 
accelerated  method,  therefore,  is  the  method  proposed  in  Qlanry  (1976),  Manry- 
Aggarwal  (1976^  ; we  call  it  the  Manry-Aggarwal  method. 


3.  Improving  the  step  size  t^ 

Although  the  Manry-Aggarwal  choice  of  the  step  size  t^  dramatically  improves 
convergence  as  compared  with  the  unaccelerated  basic  method,  further  improvement  seems 
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possible  at  little  cost.  Suppose  that  we  wish  to  move  in  a direction  g 
(1) 


(1) 


from 


b'*'  to  by  +t^  at  present  we  have 

but  in  the  next  Section  we  will  use  different  directions.  In  the  frequency  domain, 
we  will  then  have 


» ■ i» 

where  is  computed  from  as  is  from  namely  by  Equations 

1.1  or  1.4.  We  will  try  to  choose  t^  so  as  approximately  to  minimize  f,  that  is, 

(3.1)  g(t)  = E E (\^-  + tPfci  I)  . 

k“o  i“o 

as  a function  of  t.  This  function  is  difficult  to  handle  because  of  the  term 

+ tP^^^ I ; writing  this  for  the  momemt  as  |H+tP|,  we  derive  a simple  approx- 
imar ion. 

We  have  |H  + tPl  - \/  (H  + tP)  (H  + tP)  , where  " denotes  complex  conjugation. 

Therefore  |H  + tP|  = [HH  + t (HP  + HP)  + t^PP]**  - [ | H | ^ + t (HP  + HP)  + t^|P|^]^  = 

|h|(1+  t(HP+HP)/lH|^+t^|P|^/|Hfl'*.  Since  we  typically  expect  the  steps  and 

tP^^^  to  be  small,  we  can  approximate  the  square  root  (1+x)  above  by  the  terms 

" jj  2 

in  t through  the  second  power  in  the  power  series  (1+x)  “l  + l/2x-  1/8  x +... 

After  some  rearrangement,  this  gives 

(3.2)  |H  + tP|  = |h|  +t(PH  + PH)/(2|H|)  + t^ 

For  each  term  in  Equation  3.1,  by  substituting  Equation  3.2  and  rearranging  we 
therefore  have 

(A  - |H  + tP|)^  - A^-2AlH  + tPl  + |H  + tpl^ 
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lP| 

2Th1 


(PH  + PH)‘ 

8|hP 
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- A^-2A|H  + tP|  + |h|^  + t(PH+HP) 

= (A-  |H|)^  + t(PH  + PH)(l  - 


+ t^|p|^ 


(I-  ^)|P|^ 


+ A 


(PH  + PH)^ 

4|H|^ 


Substituting  this  approximation  in  Equation  3.1  gives  us  a quadratic  polynomial 
at  t,  whose  minimum  is  easily  found  to  occur  at 


2|P 


(i)|2 
k£  ' 


(3.3) 


1 1 - 

. . Vki.  ”kd  ^ \l  J 

1 l<rli 

This  expression  is  indeed  formidable,  but  It  really  is  easy  to  compute;  the  work 
is  small  compared  to  the  cost  of  the  two  discrete  Fourier  transforms  required  at 
each  iteration.  If  this  choice  of  t^  does  not  decrease  f,  then  we  instead  use 
^1  ” this  step-size,  we  require  about  75%  of  the  iterations  necessary 

with  the  Manry-Aggarwal  method  to  reduce  the  error  f to  a given  amount.  We  show 
no  numerical  examples  here,  however,  because  still  better  improvements  will  be  made 
in  the  next  section. 

4.  Improving  the  search  direction 

As  we  remarked  in  Section  2,  the  direction  used  In  the  Manry-Aggarwal  method 
is  simply  the  steepest  descent  direction  for  the  function  f of  Equation  2.4. 

Steepest  descent,  however,  is  no  longer  widely  used  by  people  in  numerical  opitmization* 
various  conjugate  gradient  or  conjugate  direction  methods  are  much  more  successful. 

The  so-called  variable-metric  versions  of  these  methods  are  unattractive  because 
they  require  storage  in  proportion  Co  Che  square  of  Che  number  In  our  case) 


of  variables.  We  therefore  opt  for  the  version  of  conjugate  gradients  described 
in  [Fletcher-Reeves  (1964),  Poljak  (1969),  Polak-Riblere  (1969)]. 

In  this  method,  the  Initial  direction  Is  the  steepest  descent  direction 

-Vf(h^^^).  Thereafter  we  take 


(4.1) 


- - f(h^^^  + 


bi_l  “ <Vf(h^^^-  Vf(h^^"^^,  Vf(h^^^>/IIVf(h^*■^^|P  . 


where  (•>•)  denotes  the  usual  inner  product 


Mj^-1  Nj^-1 


(s.x)  = E E 


k”0  l-o 


Again,  this  changed  direction  costs  little  to  compute  in  comparison  with  the  two 

discrete  Fourier  transforms  per  iteration. 

Thus  our  final  algorithm  improving  the  Manry-Aggarwal  method  consists  of  using 

the  search  direction  given  by  Equation  4.1  in  conjunction  with  the  step  size  given 

by  Equation  3.11;  whenever  f (h^^^  + t .p^^^ ) is  not  less  than  f(h^^^),  we  revert 

* - 

to  the  basic  iterative  step  and  let  the  same  as  taking  a step 

^i  * 2^  steepest  descent  direction.  This  method  dramatically  Improves  the 

Manry-Aggarwal  method,  as  we  show  in  the  next  section.  Before  considering  numerical 
experiments,  we  must  compare  the  storage  required  by  the  two  methods  now  under  con- 
sideration. The  Manry-Aggarwal  requires  on  the  order  of  2MN  + 3Mj^Nj^  locations 
(plus  terms  of  lower  order).  Because  of  needs  to  store  previous  gradients,  et  cetera, 
our  new  method  requires  on  the  order  of  4MN + 6M^^Nj  locations.  Thus  the  new  method 


requires  about  twice  the  storage. 
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5.  Some  simple  numerical  experiments 

We  show  here  a comparison  of  the  Manry-Aggarwal  method  against  our  method  on  two 

problems  from  [Manry  (1976)].  In  both  cases  we  take  M = N = 32  and  Thus 

there  are  64  unknown  coefficients  h to  determine  in  an  attempt  to  match  1024 

mn 

amplitudes  Manry-Aggarwal  method  requires  about  4500  locations  of  storage 

compared  with  about  9000  for  our  method. 

In  all  cases  we  start  the  algorithm  by  letting  o£klM-l, 

o £‘^N-1,  we  compute  h'’  from  * via  Equation  1.2,  and  finally  we  truncate 

to  let  We  show  the  results  of  100  Iterations  on  the  old  method  and 

as  many  iterations  of  the  new  method  as  needed  until  the  change  in  f in  one 
iteration  is  less  than  10  and 

Example  1 

The  desired  amplitude  response  is  a ring.  More  precisely. 


1 if  1 < 0,^  + 0^  £ 4 for  o £ k £ 31,  o £ £ 31 

k Jc 

0.2  otherwise 


where 


kit 

16 

(k-32)  TT 
16 


for  o ^k  £15 
for  16  £k  £ 31 


With  h^°^ 


determined  as  above,  initially  we  have 


“ .09575649,  ||Vf(h^°^)||  - 233. 


After  100  iterations  (200  Fourier  transforms)  the  Manry-Aggarwal  algorithm  has 
reduced  f to  a value  of  .024544519  and  ||Vf||  to  approximately  .38.  The  new 
algorithm  accomplishes  this  for  f in  52  transforms,  at  which  point  f is  reduced 
to  .024541657  and  ||  Vf  ||  to  .48;  after  78  transforms  we  have  f • .02454094819  and 
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llVfil  » .0095,  while  110  transforms  yield  f - .02A54094777  and  ||Vf  II  = .000079. 
More  data  are  displayed  in  Tables  5.1  and  5.2  below.  Both  methods  reduce  f to 
a reasonable  value  in  a few  iterations,  but  the  new  method  improves  the  f-values 
much  better  thereafter. 


10  20  30  40  50  60  70  80  90  100 

1/2  number  of  discrete  Fourier  transforms 


Table  5.1.  |^7f||  for  Example  1. 


10  20  30  40  50  60  70  80  90  100 

1/2  number  of  discrete  Fourier  transforms 


Table  5.2  f for  Example  2. 
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Example  2. 

The  desired  amplitude  response  Is  a diamond.  More  precisely. 


1 if  |0kl  + |ej  2 

o otherwise 


where  9 Is  as  In  Example  1.  With  determined  as  above,  initially  we  have 


f(h^°^  = .08179448,  ||Vf(h^°^)||  = 257. 


After  100  iterations  (200  Fourier  transforms)  the  Manry-Aggarwal  algorithm  has 
reduced  f to  a value  of  .0186557567  and  ||Vf  ||  to  approximately  3.9.  The  new 
algorithm  accomplishes  this  for  f in  52  transforms  at  which  point  f is  reduced 
to  .018574098  and  ||Vf|l  to  14;  after  120  transforms  we  have  f = .01800420801  and 
||Vf||  « .01,  while  156  transforms  yield  f = .01800420785  and  ||Vf||  = .00007.  More 
data  are  displayed  in  Tables  5.3  and  5.4  below. 


o 

X 


X “ New  method 
o » Manry-Aggarwal  method 
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018  10 


.-4 


10 


r5 


10 


.-6 


10 


20  30  40  50  60  70  80  90  100 


1/2  number  of  discrete  Fourier  transforms 


Table  5.4  f for  Example  2 


6.  References 


1.  Daniel,  James  W.  (1971),  The  approximate  minimization  of  functionals,  Prentice- 
Hall,  Englewood  Cliffs,  New  Jersey. 

2.  Fletcher,  R. , C.  Reeves  (1964),  "Function  minimization  by  conjugate  gradients," 
Comput.  J. , vol.  7,  149-154. 

3.  Manry,  Michael  T.  (1976),  "The  design  of  multidimensional  FIR  digital  filters 
by  phase  correction,  dissertation  for  the  Ph.D.  in  Electrical  Engineering, 

The  University  of  Texas  at  Austin,  Austin,  Texas. 

4.  Manry,  Michael  T.,  J.K.  Aggarwal  (1976),  "The  design  of  multidimensional  FIR 
digital  filters  by  phase  correction,"  IEEE  Trans.  Circuits  and  Systems,  vol. 
CAS-23,  185-199. 

5.  Polak,  E.,  G.  Ribierc  (1969),  "Note  sur  la  convergence  de  methodes  de  directions 
conjugees,"  Rev.  Fr.  Inform.  Rech.  Open.  vol.  16-Rl,  35-43. 

6.  Poljak,  B.T.  (1969),  "Conjugate  gradient  method  in  extremal  problems,"  J.  Comp. 
Math,  and  Math.  Phys.,  vol.  9,  809-821  (Russian). 


un.cla£al£le.d  

Vf  mifv  rlBSsificnfion 


OOCUMiHl  (OMROL  DAIA  -R«iD 

of  Korfv  of  oftAfforf  O'vo  • m*notmUttn  m-iml  bm 


t O'ViriN  « t IN  n ACVlVIty  o«trf»or^ 


The  University  of  Texas  at  Austin 


>•  Rorumiv  r i *««io«rAi«o*< 

unclassif ied  _ 


Improving  the  Manr y-Aggarwa 1 Method  for  Designing  Multi- 
Dimensional  FIR  Digital  Filters  ' 


R n0^opif*tlV^  (1%  f*»  of  mrtf4  vo  r*«im  ) 

Center  for  Numerical  Analysis 


ames  W. /Daniel 


N0OO14-67-A 


A V A a AMILIt  Y ' UtMir  A T(OH  NOYfCCt 


I AA  OTMCM  aiO<Y)  fAr««-  *o«v  ft#  •« 

rAfrorO 


Approved  for  public  release;  distribution  unlimited 


l>  4lirP(.  tMtNTAItV  NOTt« 


It  (PONSOniNO  MILITABY  ACTIVITY 


Mathematics  Branch,  Office  of 
Naval  Research,  Washington,  D.C. 


It  ASITBACT 

^ M.T.  Manry  and  J.K.  Ag^rwal  recently  described  an  algorithm  for 
use  in  the  design  of  mul ^1-d Imens lonal  FIR  digital  filters  by  phase 
correction.  As  they  observe,  their  method  can  be  viewed  as  the 
steepest  descent  method  for  minimizing  a certain  function  f (x) : 


given  an  approximate  solution  x'^y  a new  approximation  is  “ 

t^J^^  where  chosen  by  a simple  rule.  We 

derive  here  an  Improved  rule  for  determining  t^  and  an  Improved 
direction  (essentially  the  Fletcher-Reeve  s^^on  j u gat  e-grad  lent 

direction).  The  resulting  method  appears  to  be  two  to  three  times 
as  fast  as  the  Manry-Aggarwal  method;  the  additional  cost  is  pri- 
marily in  storage,  which  roughly  doubles.  •/ 

Key  Words:  FIR  digital  filters;  filter  design;  ^hase  correction; 

Manry-Aggarwal  method. 
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